This tutorial … Scottish Remote Sensing Portal.
library(lidR)# process LAZ files
library(raster)# deal with raster
library(rayshader) # 3D viz
library(rgl) # interactive plot
library(sf) # spatial classmotte_point = st_read("data/canmore_point.shp", quiet = TRUE)
motte_buffer = st_buffer(motte_point,dist = 50)Next, I am reading the LAZ files and clipping the point cloud to the extent of immidiate area around the motte.
#read laz files
las = readLAS("data/NX6055_4PPM_LAS_PHASE3.laz")
motte = clip_roi(las, motte_buffer)
#plot lidar point cloud
plot(motte, bg = "white")
rglwidget()grab and rotate the model !
In this step I am running standard alghoritms from LiDR package to compute Digital Surface Model (DSM) and Digitial Terrain Model (DTM).
rgl.clear()
# create dsm and dtm
dsm = grid_canopy(motte, 0.5, pitfree())
# assign coordinate system
crs(dsm) = CRS("+init=epsg:27700")## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum OSGB_1936 in CRS definition
st_crs(27700)## Coordinate Reference System:
## User input: EPSG:27700
## wkt:
## PROJCRS["OSGB 1936 / British National Grid",
## BASEGEOGCRS["OSGB 1936",
## DATUM["OSGB 1936",
## ELLIPSOID["Airy 1830",6377563.396,299.3249646,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4277]],
## CONVERSION["British National Grid",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",49,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",-2,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",0.9996012717,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",400000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",-100000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["unknown"],
## AREA["UK - Britain and UKCS 49°46'N to 61°01'N, 7°33'W to 3°33'E"],
## BBOX[49.75,-9.2,61.14,2.88]],
## ID["EPSG",27700]]
# create dtm
dtm = grid_terrain(motte, 0.5, algorithm = knnidw(k = 6L, p = 2))
# addign coordinate system
crs(dtm) = CRS("+init=epsg:27700")## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum OSGB_1936 in CRS definition
par(mfrow = c(1,2))
plot(dtm, main = "DTM", col = height.colors(50))## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Airy 1830 ellipsoid in CRS
## definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Airy 1830 ellipsoid in CRS
## definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Airy 1830 ellipsoid in CRS
## definition
plot(dsm, main = "DSM",col = height.colors(50))## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Airy 1830 ellipsoid in CRS
## definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Airy 1830 ellipsoid in CRS
## definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Airy 1830 ellipsoid in CRS
## definition
# dsm
slope_dsm = terrain(dsm, opt = 'slope')## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Airy 1830 ellipsoid in CRS
## definition
aspect_dsm = terrain(dsm, opt = 'aspect')## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Airy 1830 ellipsoid in CRS
## definition
hill_dsm = hillShade(slope_dsm, aspect_dsm, angle = 40, direction = 270)## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Airy 1830 ellipsoid in CRS
## definition
# dtm
slope_dtm = terrain(dtm, opt = 'slope')## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Airy 1830 ellipsoid in CRS
## definition
aspect_dtm = terrain(dtm, opt = 'aspect')## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Airy 1830 ellipsoid in CRS
## definition
hill_dtm = hillShade(slope_dtm, aspect_dtm, angle = 5, direction = 270)## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Airy 1830 ellipsoid in CRS
## definition
#plot
par(mfrow = c(1,2))
plot(hill_dtm, main = "DTM Hilshade", col = grey.colors(100, start = 0, end = 1),
legend = FALSE)## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Airy 1830 ellipsoid in CRS
## definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Airy 1830 ellipsoid in CRS
## definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Airy 1830 ellipsoid in CRS
## definition
plot(hill_dsm, main = "DSM Hillshade", col = grey.colors(100, start = 0, end = 1))## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Airy 1830 ellipsoid in CRS
## definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Airy 1830 ellipsoid in CRS
## definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Airy 1830 ellipsoid in CRS
## definition
#And convert it to a matrix:
elmat = raster_to_matrix(dtm)
elmat %>%
sphere_shade(texture = "imhof1") %>%
add_shadow(ray_shade(elmat, zscale = 0.5, sunaltitude = 30,lambert = TRUE),
max_darken = 1) %>%
add_shadow(lamb_shade(elmat,zscale = 0.5,sunaltitude = 30), max_darken = 0.2) %>%
add_shadow(ambient_shade(elmat, zscale = 0.5), max_darken = 0.2) %>%
plot_3d(elmat, zscale = 0.5,windowsize = c(1000,1000),
phi = 40, theta = 180, zoom = 0.8, fov = 1)
rglwidget()